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Abstract 

We present test results for the smooth lattice method using an Oppenheimer- 
Snyder spacetime. The results are in excellent agreement with theory 
and numerical results from other authors. 



1 Introduction 

In recent times many numerical relativists have good reason to celebrate - the 
long battle to secure the holy grail pQ is over (though some might prefer to 
redraw the battle lines) . The works of Pretorius [21 E] and others (H [5] have 
opened a new era for computational general relativity. This has spawned 
many new projects that directly address the needs of the gravitational wave 
community. Many groups are now running detailed simulations of binary 
systems in full general relativity as a matter of course. Does this mean that 
the development of computational methods for general relativity is now over? 
The experience in other fields would suggest otherwise, look for example at 
computational fluid dynamics where a multitude of techniques are commonly 
used, including spectral methods, finite element methods, smooth particle 
hydrodynamics, high resolution shock capture methods and the list goes 
on. The important point to note is that one method does not solve all the 
problems and thus in numerical relativity it is wise, even in the face of the 
current successes, to seek other methods to solve the Einstein equations. 
It is in that spirit that we have been developing what we call the smooth 
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lattice method [6J [7J [8] . This is a fundamentally discrete approach to general 
relativity based on a large collection of short geodesic segments connected 
to form a lattice representation of spacetime. The Einstein equations are 
cast as evolution equations for the leg-lengths with the Riemann and energy- 
momentum tensors acting as sources. Of course the Riemann tensor must be 
computed from the leg-lengths and this can be done in a number of related 
ways, such as by fitting a local Riemann normal coordinate expansion to a 
local cluster of legs or to use the geodesic deviation equation, or, and with 
more generality, to use the second variation of arc-length. Past applications 
of the method have included a full 3+1 simulation of the vacuum Kasner 
cosmology [8] and a 1+1 maximally sliced Schwarzschild spacetime [6]. In 
both cases the simulations were stable and showed excellent agreement with 
the known solutions while showing no signs of instabilities (the maximally 
sliced Schwarzschild solution ran for t > 1000m and was stopped only because 
there was no point in running the code any longer). 

In this paper we report on our recent work using the Oppenheimer-Snyder 
[9] spacetime as a benchmark for our smooth lattice method [6J [8] . We 
chose this spacetime for many reasons, it has been cited by many authors 
[TUl [TTJ [T2J EH El H5J US] as a standard benchmark for numerical codes (and 
thus comparative results are available), the analytic solution is known (in a 
number of time slicings), the equations are simple and there are many simple 
diagnostics that can be used to check the accuracy of the results (as described 
in sections 



11. 12) 



In an impressive series of papers, Shapiro and Teukolsky ( [TUl EH EH]) used 
the Oppenheimer-Snyder spacetime as the first in a series of test cases. They 
were motivated by certain problems in relativistic stellar dynamics (such as 
the formation of neutron stars and black holes from supernova) and they 
developed a set of codes based on the standard ADM equations, adapted to 
spherical symmetry, in both maximal and polar slicing and using an iV-body 
particle simulation for the hydrodynamics. They made limited use of the 
exact Schwarzschild solution to develop an outer boundary condition for the 
lapse function while using both the Schwarzschild and FRW solutions to set 
the initial data. Though their discussion on the size of their errors is brief 
(for the Oppenheimer-Snyder test case), they did note that the errors were 
of the order of a percent or so (for a system with 240 grid points and 1180 
dust particles). In a later work, Baumgarte et al. [15] extended their work 
by expressing the metric and the equations in terms of an out-going null 
coordinate. This leads to a slicing that covers all of the spacetime outside 
(and arbitrarily close to) the event horizon. In this version of their code 
Baumgarte et al. [15] chose to solve only the equations for the dust ball by 
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using the Schwarzschild solution as an outer boundary condition. 



This idea, to replace the exterior equations with the known Schwarzschild 
solution, has been used by Gourgoulhon [13], Schinder et al. [12] and Romero 
et al. [16]. Gourgoulhon [13j used a radial gauge and polar slicing while 
solving the equations using a spectral method and reported errors in the 
metric variables between 10~ 7 to 10~ 5 . However, with the onset of the Gibbs 
phenomena, the code could only be run until the central lapse collapsed to 
around 2 x 1CT 3 . Schinder et al. [12] used the same equations as Gourgoulhon 
[13] but with a discretisation based on a standard finite difference scheme. 
They reported errors of order 1% for evolution times similar to those of 
Gourgoulhon [13] . The work of Romero et al. [16] differs from that of Schinder 
et al. [12] in that they used high-resolution shock capture methods for the 
hydrodynamics. They report evolutions down to a central lapse of 1.3 x 10~ 10 . 

Our results compare very well against those given above with our errors 
being of the order fractions of a per-cent for 1200 grid points. Our code 
runs, without any signs of instabilities, for maximal slicing out to t = 500m 



where the central lapse has collapsed to 10~ 110 (see Figure 27). We make no 
use of the known solutions other than the conservation of local rest energy 
(we use a particle like method to compute the rest density). We also provide 



extensive comparisons of our results with the exact solution (see section 12). 



In the following sections we will describe all aspects of our code, including 
the design of the lattice (section [2]), the curvature and evolution equations 
(sections |3j [4] and |6| , computing the density (section [7]), the junction condi- 



tions (section [9]), setting the initial data (section 10.5) and finally the results 
(section [l2]). 



We will make frequent reference to two papers, our earlier paper on the 
Schwarzschild spacetime [8] and a companion paper showing how the Einstein 



equation can be applied to the lattice [17]. We will refer to these as Paper 1 



and Paper 2 respectively 



2 The Oppenheimer-Snyder lattice 

What design should we choose for the lattice? We will take a minimalist 
approach - build the simplest lattice that captures the required symmetries 
while being sufficiently general to allow the full dynamics to be expressed 
through the evolution of the lattice data. Here is a construction of such 
a lattice. Take a single spacelike radial geodesic, in one Cauchy surface, 
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extending from the centre of the dust ball out to the distant asymptotically 
flat regions and sub-divide it into a series of short legs with lengths denoted 
by L zz . We will refer to the end points of each leg as the lattice nodes. Note 
that we are free to choose the L zz as we see fit (in the same way that we are 
free to choose the lapse function in an ADM evolution). Now construct a 
clone of this geodesic by rotating it through any small angle (while remaining 
in the Cauchy surface). Finally connect the corresponding nodes of the pair 
of geodesies by a second set of geodesic legs, with lengths denoted this time by 
L xx (see Figure [TJ. We now have a spacelike 3-dimensional lattice contained 
within one Cauchy surface. From here on in we allow this lattice to vary 
smoothly with time. 

Note that each leg of this lattice is a geodesic segment of the 3-metric of the 
Cauchy surface. We could also connect the nodes of the lattice with geodesic 
segments of the full 4-dimensional spacetime (much like constructing chords 
to arcs of a circle). This gives us two representations of the lattice, both 
sharing the same node points with the first composed of short 3-geodesics 
and the second composed of short 4-geodesics. Suppose that typical leg- 
lengths in the two representations are 3 Lij and 4 Ljj respectively. Then it is 
not hard to see that 4 Lij = 3 Lij + 0( 3 L^). The upshot is that in all of our 
equation in this paper we are free to use either representation (the differences 
being at least as small as the truncation errors). 

The L zz and L xx are all that we need to describe the geometry of each Cauchy 
surface but we also need some way to represent the dust ball on the lattice. 
Again, we shall take a minimalist approach - we know that the dust can be 
described as a set of particles travelling on timelike geodesies with conserved 
rest mass. Thus we add a series of dust particles on the radial geodesic with 
each particle carrying a conserved rest mass. 

As noted above, we are free to distribute the lattice nodes as we see fit. How 
should we do this? We know that the dust ball will collapse so it makes 
sense to tie the lattice nodes to the dust particles, i.e. the lattice nodes 
follow in-falling timelike geodesies. But what of the nodes outside the dust 
ball? Again, by appeal to simplicity, we demand that every lattice node, 
interior and exterior, follow the in-falling timelike geodesies. In this scheme 
the lattice nodes do not follow the trajectories normal to the Cauchy surface 



in contrast to the scheme in Paper 1 ). This introduces a drift vector 7^ (see 



Figure [3]) (which is similar to but distinct from the shift vector, see j8]). 



The lattice just described differs from the Schwarzschild lattice of |Paper~T 



in a number of important ways - it contains an internal boundary (the edge 
of the dust ball), the lattice nodes are not at rest in each Cauchy surface, 
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the lattice carries a set of dust particles and at the inner boundary L xx = 0. 
Thus we will need to develop new boundary conditions (section [9]), new 
evolution equations for the nodes (i.e. adapt the geodesic equations to the 
lattice, section [6]) and an algorithm to compute the rest energy density from 
the rest masses carried by the dust particles (section [7]). 



In Paper 1 we employed Riemann normal coordinates as a stepping stone 
to develop the purely scalar equations for the leg-lengths, time derivatives, 
constraints etc. We went on to speculate whether or not these coordinates 
imbued the numerical scheme with any favourable properties (we argued that 
they did not). One way to avoid this coordinate issue is simply to derive the 
equations without reference to a coordinate system. In this paper we will 
represent tensors, such as the Riemann and extrinsic curvatures, by their 
frame components. We will use an orthonormal frame built as follows. We 
choose the first two basis vectors, and m^, to be the unit tangent vectors 
to L xx and L zz respectively at the mid-point of L xx , see Figure ([!]). The 
remaining two basis vectors (m^ and m^) can be chosen freely (subject to 
the orthornormal condition, e.g. could be chosen as the unit normal to 
the Cauchy surface). With this choice of basis a typical frame component for 
the extrinsic curvature could be written as K^m^m^. Such notation quickly 
becomes tiresome so we will introduce the abbreviation JC a b to represent 
Kfa/fTi^m^ with an obvious generalisation to other tensors. 

We will allow a slight variation to this notation. On occasions we will find it 
useful to refer to a leg by its end points, such as i and j. That leg will have 
its own unit tangent vector which we denote by m M . We will then take K-ij 
to be K tlu m fJ- m u . This small change will only ever be used for the extrinsic 
curvature. 

The dust particles follow future pointing timelike geodesies. We will use 
to denote the velocity 4-vector of the dust particles and we will record the 
frame components as v n = —v^n^ and v z = v^m^, where is the unit normal 
to the Cauchy surface. 



The notation just introduced sits quite nicely with the notation used in |Paper 
[TJ In that paper we wrote K xx , for example, to denote the x — x coordinate 
components of K^ v in the local Riemann normal frame. In that frame we 
chose the three metric g^x) at the origin to equal diag(l, 1, 1) and the basis 
vectors mf^ to have values 6%. Thus K xx = K^m^m^. The upshot is that 



coordinate components of Paper 1 have the same numerical values as the 



frame components used in this paper. Thus we would reasonably expect 



that the equations used in Paper 1 should carry over to this paper with only 



minor changes to accommodate the introduction of the dust. This indeed 
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proves to be the case (which is reassuring). The details will be presented in 
section [4] where we will use a formalism developed in Paper 2 to derive, from 
scratch, the evolution equations for the lattice. 



3 The Riemann curvatures 



The question here is: How do we compute the Riemann curvatures, lZ xyxy and 
TZ-xzxz, from the L xx and L zz l In our previous paper, |Paper 1[ we computed 
two Riemann curvatures, TZ xyxy and 7Z XZXZ , using 

= — — — h TZ XZXZ L XX geodesic deviation (3.1) 



dz 2 

= d(L 2 xx n xyxy ) _ Uxzxz d^L Bianchi identity (32) 

dz dz 

These equations will be used as follows. First we use the geodesic deviation 
equation to compute the 1Z XZXZ for each node across the lattice. This then 
allows us to integrate the Bianchi identity for lZ xyxy from the centre to the 
outer boundary. 

This scheme sounds simple but there a number of (obvious) complications. 
Firstly, the equations are singular at the centre (where L xx = 0) and secondly, 
1Z XZXZ will not be continuous across the junction. These complications are 
new to this investigation but we also inherit one further complication from 
|Paper 1 what boundary condition should we use at z = when integrating 



the Bianchi identity? This last problem is rather easy to deal with. At the 
centre of the dust ball we know that the metric must be isotropic and thus 
we can be certain that lZ xyxy = 1Z XZXZ at z — 0. 

How do we handle the singularity at z = 0? Again, by symmetry arguments 
we can assert that 1Z XZXZ ^ = at z = 0. Thus in a small neighbourhood 
of z = we must have TZ XZXZ = A + Bz 2 where A and B are independent 
of position and z is the radial proper distance measured from z — 0. Thus 
it is not unreasonable to use a quadratic interpolation of 1Z XZXZ to estimate 
T^-xzxz at z = 0. Our experience shows that this works very well but it does 



require some care (see section 10.1 for the full details). Before dealing with 
the junction issue we should emphasise that this process is an interpolation 
rather than an extrapolation of the data to z = 0. To see this just imagine 
extending the radial geodesies of the lattice through z = so that we can 
use 7Z XZXZ = A + Bz 2 with z in a range — z < z < zq for some small z . 
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The frame components 1Z XX = R^m^m^ and 1Z ZZ = R^m^mK, of the Ricci 
tensor are rather easy to construct from lZ xyxy and 1Z XZXZ . Using the or- 
thonormal frame m%, and we can easily deduce that 

T^xx T^xyxy ~\~ T^xzxz (^"^) 
^"(/y T^xyxy ~\~ T^-xzxz (3'4) 
T^zz T^xzxz (3-5) 

with all other lZ a b = 0. From these equations it is easy to verify that the 
scalar curvature is given by 

R '2nRs X y X y -\- 4:'TZ. XZXZ (3.6) 

The one remaining complication is the discontinuity in TZ XZXZ at the junction. 
This will be discussed in detail in section |9l 



4 The evolution equations 



The equations of Paper 1 have served us well so far but now we must chart 
a new path. The reason is that, unlike our approach in Paper 1 here we 



allow the lattice nodes to drift across the Cauchy surfaces and this will in- 
troduce extra terms in the evolution equations. There is also the issue of 
introducing the energy momentum sources but, as we shall see later, this is 
really very easy to do (it amounts to little more than adding a term of the 
form SnkT^m^m 1 ' to the vacuum equations). So how do we develop evolu- 



tion equations for a non-zero drift vector? In Paper 2 we showed how the 
standard 3+1 ADM equations with a zero shift vector can be recovered from 
the equations for the second variation of arc-length. And lengths of 

geodesies are central to our smooth lattice approach this new formalism is 
well suited to our current task. 



We begin by recalling from |Paper 2 the equations for the first and second 
variation of the geodesic segment that connects nodes i and j 



dL 



dt 



m (1 m l 'f. l/ ds 



(4.1) 



+ 



(4.2) 
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It is tempting to jump in by setting = Nri^ + 7^ and to let the equations 
take us where they will. Indeed this works well for the first variation. We 
start be making the said substitution and massage the result as follows 

= K(iVn„)K + KtJ5 

j 

m^m u {Nn^). p ds + [m M 7 (U ]^ 

' m»m u (-NK^) ds + [m^l 
= -NKyLy + KtmU 

In the second last line we have used Nn^. v = — (_LiV M ) n u — NK^ U and 
m^n^ = while in writing the last line we have assumed that the leg-length 
is sufficiently short that the integrand can estimated by a simple quadrature, 



in this case the mid-point rule (later in section 10.3 we will have reason to 
change this to a Trapezoidal rule). A similar equation can be found in Paper 
[2] (differing only in the absence of the 7 terms). For the pair of legs L xx and 
L zz we thus obtain 

= -NK, XX L XX + [m^] 

^ = -NK ZZ L ZZ + K 7 ,l 

and to keep the notation a little less cluttered we have not written the end 
points on the [• • • ] terms. 

What can we say about the [m fl , y fl ] terms? Let 7^ be the unit vector parallel 
to 7^. Then we can immediately use the first variation equation once again 
(see Figure ^ to deduce that [m 2 . M 7 M ] equals dL xx /dz where, as usual, z is 
the radial proper distance measured along G\. However, 7 M = 7 M / r y z and 
by spherical symmetry we know that r ) z does not change from one radial 
geodesic to the next. Thus we deduce that [m x ^] = ^ z dL xx jdz. We now 
turn to the other leg, L zz . In this case and 7^ are parallel and thus we 
can not invoke the first variation equation. But that is of no concern simply 
because m z ^ = 7 2 = Nv z /v n . Thus we have [m z ^y z ] = [Nv z /v n ]. The 
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equations for the first time derivatives of L xx and L zz can now be written as 



dL n 



dt 
dL z , 
dt 



-NJC ZZ L ZZ + 



Nv,\ dL 



dz 



(4.3) 
(4.4) 



We turn our attention now to adapting the equations for the second variation 
to our simple lattice. 

Our job would be greatly simplified if it happened that 7 M = but on the 
current lattice that is not the case. So we now introduce a second lattice 
on which we set 7^ = 0. The nodes of the first lattice will follow the dust 
particles while those of the second lattice will follow trajectories normal to 
their Cauchy surfaces (which will differ from those of the first lattice). Note 
that the second lattice has been introduced solely to aid the exposition - the 
second lattice will never be needed nor used in our actual computer programs. 
To keep the bookkeeping clear we will identify data on the second lattice by 
the addition of a dash. The second lattice is created at some generic time, 
say t = to, and we choose to assign identical initial data to both lattices, i.e. 



L\ - = Lij, K'^ = /Cjj etc. on S(t )- We have no reason to use distinct Cauchy 
surfaces for each lattice (we only want to set 7^ = 0) so we are free to set 



N' = N and dN' /dt = dN/dt. It follows that we also have A R' 



across £(£q)- Our task now is to adapt the equations for the first and second 



variations to the second lattice. This has already been done in Paper 2 where 
we have shown that 



dU 



1 - -NK'^mWL^ 



dt 

d 2 L' tJ _ 1 dNdL'y 



dt 2 



N dt dt 



dJ' \ 2 

-£) +N*Kl a K^m«mPL> 



+ NN ;a/3 m a m l3 L' ij - N 2 ( 4 R^) 



m>*m u n a n (3 LL 



Where we head next depends upon what type of equations we wish to work 
with. We can develop either a second-order set of equations involving both 
dL'^/dt and d 2 L' ij /dt 2 or a first order system involving dL[-/dt and dJC'^/dt. 
We will take the second approach for two reasons, it mimics the standard 
ADM approach and, more importantly, it eliminates the dN/dt term (which 
would add undue complexity when using maximal slicing) . Between this pair 
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of equations we can easily eliminate d 2 L ; i j/dt 2 with the following result 



dIC' 



dt 



- ~ -N. al3 m a m p + N ( 4 R mu p) m ft m v n a n p 



2N (iC'^-NK^K'^m 



where /CL := K'm^m". This last equation controls the evolution of /C,- for 
the second lattice. Can we use this information to deduce the evolution of 
fCij on the first lattice? Yes, by simply splitting the evolution into a part 
parallel to the normal plus a part parallel to the drift vector. Since /C^ is a 
scalar function we can use a standard chain rule to write 

dt dt + 

and as JCij = fC'^ on £(io) we arrive at 

= -N, a pm a m p + N ( 4 R„ au p) m^m u n a n p 
+ 2N (K l3 f - NK m K^m a mP + JC^Y 

In a moment we will apply this equation to fC xx = K^ v m^.m v x and K zz = 
^nu m z Tn z but first we recall that both 7^ and m M are tangent to £, T^ u = 
pv^v v and for our spherically symmetric lattice, K^ v is diagonal and j z = 
Nv z /v n . We will also need the contracted Gauss equation, namely, 

We then find that the above equation for dK^jdi when applied to K xx and 
K zz leads to 

' /A " -N xx + N (TZ XX + KK XX - 4nkp) + (^) 1C XX}Z (4.5) 



dt 
dfC z 
~dT 



-N xx + N (1Z ZZ + KK ZZ + (4 - 8v 2 n ) nkp) + (^A K z 



(4.6) 



This pair of equations coupled with ( 4.3|4.4 ) are the evolution equations for 
the lattice. 



5 The constraints 



The general form of the Hamiltonian and momentum constraints are 
R + K 2 - K^K^ = IQ-KkT^n" 
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where K = K^u. It is a simple matter to apply these equations to the 



Schwarzschild spacetime, see Paper 1 for details. For the present case we need 
to account for the non-zero in the interior of the dust ball. We can easily 



adapt the equations of Paper 1 by simply adding on the terms 87ikT^ u n fl n u for 



the Hamiltonian and SirkT^n^m^. for the momentum constraints (projections 
in the other two directions and yield the trivial equation = 0). This 
leads to 



Ti-xyxy ^T^-xzxz ^xx ^^xx^zz 
1 / . „ 







dz 



d{L xx JC 2 
dz 



\ixkpv\ 
Airkpv n v z 



(5.1) 
(5.2) 



(for simplicity we have cleared a common factor of 2 from both equations) 



6 The particle equations 



Here we will derive the equations governing the evolution of the particle 
4- velocities. 

We will use the geodesic equation = v^- u v u to obtain evolution equations 
for the v n and v z components of the particle's 4-velocity. 

The computation are simple but do entail a few steps. We begin by writing 
dv n /dt as a directional derivative along t M = At> M . The Leibniz rule is then 
applied which in turn allows the geodesic equations = v^. u v u to be imposed. 
Finally, we use 

Nn KV = -±JSf tlt n v - NK^ U (6.1) 

to re-write n^ v in terms of the lapse and extrinsic curvatures. The details 
are as follows. 

dv 



1 

N 

A „ 1 



- -v z v n j^N z + jYYKnv 

But we also know that A = N/v n and = 7 Z = 7^ while j z = Nv z /v n so 
this last equation may be further reduced to just 

= + N - K » ( 6 - 2 ) 

dt v„ 



11 



The computations for dv z /dt are much the same, 
dv 

= A (t) n n" + w 2 m^) v v m zll . v 
= Xv n n^v u m ztl . u 

The term n )1 v u m Zfl . u can be computed by expanding = (m Zfl n >M ). u v l/ and 



then using (6.1) to obtain 



v. 



= m zu: ^/ + -^N z - K, 



,m l tm u 



which, when substituted into the previous equation for dv z /dt, leads to 

^ = ~v n N z + Nv z K zz (6.3) 
dt 

One simple check we can immediately apply to our equations is to ask: do 
they preserve the unit normalisation of v M ? Since we have chosen n M and 
to be unit vectors the question reduces to asking if d{—v\ + v 2 z )/dt vanishes 
for all t. From the above equations this is easily seen to be so. 



7 The density 

There are at least two ways to compute the density, either by solving the 
Hamiltonian constraint or by integrating the equations of motion for the 
dust, namely, = (pv^v u ). u . 

Recall that the Hamiltonian constraint is given by 

Tl xyxv + 2K XZXZ + K} xx + 2)C XX )C ZZ = Znkpvl (7.1) 

This equation is trivial to solve for p since on each Cauchy surface all of the 
other quantities are known. Notice that v^ — l + v^ and thus v n > 1 > 0. 

Using the Hamiltonian constraint is one of many tricks used in numerical 
relativity to coerce better stability properties from the evolution equations. 
The merits of doing so have been debated over the years and is not something 
we will delve into here. However as we are trying to establish the limitations 
of the smooth lattice method it makes sense to explore other methods to com- 
pute the density. So for our second method we turn to the energy-momentum 
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equations. From = (pv^v u ). v we learn two things (i) the dust particles fol- 
low time like geodesies, = v fJ, - u v u and (ii) the rest mass is conserved along 
the worldtube generated by the dust particles = d/dt J pdV where d/dt 
is the time derivative following the dust and dV is the proper volume in the 
dust's rest frame. We will need both equations to compute the density. 

Recall that we have chosen to tie the dust particles to the nodes of the 
lattice. As the nodes drift relative to the Cauchy surface there will be a 
non-zero boost between the rest frame of the dust and that of the Cauchy 
surface. Thus, in terms of the volume element dV on the Cauchy surface we 
have 

/ pv n dV = pv n dV 

J Co J C\ 

where C , C\ denote the intersections of a dust worldtube with a pair of 
Cauchy surfaces, one at time to an d another at a later time t\. 

The question which arises now is: how do we construct the three dimensional 
cross-sections Co from the 2-dimensional lattice? The solution is depicted in 
Figure [2] where we have simply taken the original lattice and rotated it by 
7r/2 about the central geodesic G\. This creates C and C\ as truncated 
pyramids with a square cross-section. In each of these we take the density 
to be constant. The volume of Co and C\ can be computed by elementary 
Euclidean geometry (the dust is minimally coupled to the geometry and thus 
curvature corrections can be ignored). This leads to 

{V\ = \{L ZZ \ ((Ly 4 + (L xx )i(L xx ) i+ i + (LL) m ) 

where (L xx )i and (L xx )i + i are the values of L xx at nodes % and i + 1 respec- 
tively. The previous conservation equation can now be re-written as 

3m; = (pv n ) t (L zz )i . + (L xx )i(L xx ) i+ i + (LL) m ) (7.2) 

where m 8 is the conserved rest mass along the worldtube (m, is set as part 
of the initial conditions). The v n are estimated at the centre of each cell 
by quadratic interpolation from the neighbouring nodes (which will draw in 
nodes beyond this basic cell). This equation can then be solved for p. We 
assign that p to the centre of the cell and then use quadratic interpolation 
to estimate p at the lattice nodes. 
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8 Maximal slicing 



A maximally sliced spacetime is defined to be a spacetime for which K — 
everywhere. Such spacetimes are often constructed by first setting K = 
on an initial Cauchy surface (e.g. on a time symmetric initial slice) and then 
demanding that dKj dt = throughout the evolution. For our lattice we have 



K = 2JC XX +JC ZZ and thus from the equations ( 4.5[4.6 ) we see that dK/dt = 
provided 

= 2N XX + N )XX — N (R — Ank (l + 2v 2 n ) p) 



But in |Paper 1| we showed that under spherical symmetry 

1 dL xx . w 

N xx = t—N z 

which allows us to re-write the previous equation as 

2 dL 

= N zz + -f^N z -N(R- Ank (l + 2v 2 n ) p) (8.1) 

We treat this as an ordinary differential equation for N. The boundary 
conditions are simple, at z = we require dN/dz = while at the outer 
boundary we require 1 = lim z ^ 00 N. Note also that the differential equation 
is singular at z = (due to the l/L xx term). We deal with this by appealing 
to the spherical symmetry of the solution at z = to deduce that N xx = N iXX 
and thus our original differential equation for N can be re-written as 

= 3N iZZ -N(R- Ank (l + 2v 2 n ) p) at z = (8.2) 

which is clearly non-singular. The same result can also be obtained by ap- 
plying l'Hopital's rule to (1/L XX ) (dN/dz) as z — > 0. At the junction we know 
that p and R suffer a jump discontinuity. Thus we expect a corresponding 
jump discontinuity in d 2 N/dz 2 which in turn forces both N and dN/dz to be 
continuous across the junction. This adds extra constraints to the numerical 
solution of the above equation. We will cover this in more detail in section 



10.2| but for the moment we note that our method computes two separate 
solutions, one for either side of the junction, which are then matched at the 
junction. 



9 The junction conditions 

Darmois [Tj5] and later Israel [2D] developed a very elegant approach to handle 
discontinuities in a metric in General Relativity. However, their method 
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requires some work to push through so we defer the details to Appendix [A] 
preferring instead to present here a direct approach. 



By integrating the geodesic deviation equation (3.1) over a short interval 
z G (— e, +e) we obtain 



dz 



I T^xzxzLxx dz 



If we require 1Z XZXZ to be bounded on each Cauchy surface then we must have 







lim 



dL XQ 
dz 



+ c 



and thus dL xx /dz is continuous everywhere on the lattice and, most im- 
portantly, across the junction. We also know that L xx and dL xx /dt must 



be continuous and thus from the evolution equation (4.3) we see that = 
lim^o^CxxJ-e- From here on we shall dispense with the limits on the square 
brackets and take [. . . ] to mean lim e _+o[- ■ ■ \_\. 



Applying a similar integration to the Bianchi identity leads to 
= [L 2 xx TZ xyxy \ - lim j 7Z, 



dT 2 

^^dz 
dz 



and thus 



o = [n 



xyxyj 



since L xx must be continuous every where on the lattice. 



(9.1) 

Thus we con- 



clude that lZ xyxy is continuous on the lattice. However, by inspection of the 
Hamiltonian constraint ( |5.1[ ), we see that the same can not be said for 1Z XZXZ . 
Since we know that = [L xx ] and = [K xx \ we see that continuity of the 
Hamiltonian requires 



^-■xx^zzl 



(9.2) 



We also need suitable junction conditions for the lapse function when us- 
ing maximal slicing. First we demand that the clocks of a pair of observers 
travelling close to but on opposing sides of the junction should remain syn- 
chronised throughout their journey. Thus we find that the lapse is continuous 
across the junction, = [TV]. For the first derivative we follow the method 



outlined above. Integrating the maximal slicing equation (8.1) over the short 
interval z 6 (— e, +e) leads to 



= [N : 



+ lim 

e-*0 



dL xx 



dz 



N g -N(R- Airk (1 + 2v 2 n ) p) j dz 
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and as we except all terms in the integral to be bounded (at worst) and 
L xx > we see that this requires 



= [N, 



(9.3) 



Equations (9.1), (9.2) and (9.3) constitute the full set of junction conditions 
for our lattice. Other conditions such as = [L xx ] and = [N] are trivially 
implemented in the numerical code (they require no special care). However 
we have no freedom in our data to guarantee = [dL xx /dz\. The reason 
is that all of the L xx leg lengths are subject to the evolution equations and 
we have to live with what they dictate. Of course we expect the jump in 
dL xx /dz to be small and to vanish as the lattice is progressively refined. 



10 Numerical methods 



To obtain numerical solutions of our equations we turn once again to the tech- 
niques developed in |Paper 1 We use second order accurate finite differences 
(on a non-uniform grid) for all of the spatial derivatives, such as dL xx /dz and 
d 2 N/dz 2 (though with a two exceptions, as noted below in section 10.4, for 



the the three nodes centred on the junction). The time integration employs 
a standard 4th-order Runge-Kutta method and the time step is chosen so 
that the Courant factor for the smallest L zz on the lattice is 1/2 (the leg on 
which this occurs lies on the surface of the dust ball). 

The lattice and its attendant equations in this paper differ most notably from 



those of Paper 1 by the presence of the dust ball. This not only introduces 
new terms in the equations but it also forces many of the variables, or their 
derivatives, to be discontinuous at the junction. Dealing with these discon- 



tinuities requires some care. For the geodesic deviation equation (3.1), the 



Bianchi identity (3.2) and the maximal lapse equation 8.1 the general ap 



proach is to solve those equations twice, once on either side of the junction, 
and then use the junction conditions to match the solutions. The details are 
as follows. 



10.1 The Riemann curvatures 



The discretised forms of the geodesic deviation equation (3.1 ) and the Bianchi 



identity (3.2) were given in Paper 1 and, apart form some minor notational 
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changes, are equivalent to the following pair of equations 



2 {pxx)i (R-xyxy)i = {^xx)i {(T^xzxz)i + (T^xzxz)i-l) 

"I - (-^a!x)i_i (^(J^xyxy)i—1 ^J^-xzxz)i ^J^-xzxz)i—\j 

where the second derivatives of L xx are computed using the second order 



non-uniform finite differences (as described in section 10.4). 



Our plan is to use this pair of equations to calculate the Riemann curvatures 
on the lattice but we immediately encounter two problems, the equations are 
singular at z = and, as previously noted, the second derivatives of L xx are 
not continuous across the junction. The first problem is rather easy to deal 
with. We draw upon the required spherical symmetry at z = to deduce that 
dTZxzxz/dz = at z = and thus lZ xzxz (z) = A + Bz 2 + O (z 3 ) near z = 0. 
The coefficients A and B are obtained by fitting lZ xzxz (z) = A + Bz 2 to two 
samples for (TZ xzxz )i (typically {Tl xzxz )± and (TZ XZXZ ) S for rij = 120) and then 
setting (7Z xzxz )i = A + Bz 2 for each node near z = (i.e. at z = 0,zi,z 2 
and z 3 ). For TZ XZXZ we again call on the spherical symmetry to assert that 
(R. xyX y)o = (R-xzxz)o- Our numerical experiments show that we have no need 
to use the quadratic interpolation scheme for (TZ xyxy ) near z — 0. 

We turn now to the issue of the junction. As with the lapse function, we 
compute both Riemann curvatures separately on each side of the junction. 
We first use the above equations to compute the curvatures for all of the 
interior lattice nodes excluding the node at the junction. At the junction we 
apply a series of interpolations in conjunction with the boundary conditions 
to set the curvatures on the junction and one node point outside it. The 
details are as follows. 

First we use cubic extrapolation to compute the one-sided limits lim 2 | 2j 7Z xyxy 
and lim 2 | 2j 7Z XZXZ , which we abbreviate as (H xyxy )~ and (1Z XZXZ )~. We then 



use the junction condition (9.1) and the Hamiltonian constraint (5.1), which 
we re-write as 

(J^xyxy) (7^-xyxy 

)" (10.3) 

(J^-xzxz) 2 (j^xyxy K-'xx ^xx^zz) (10.4) 

to step across the junction (with the + super-script denoting the right hand 



one-sided limit). We then return to the above discrete equations (10.1) and 
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(10.2) to compute the curvatures in the exterior region. This too requires 
some explanation. We first compute (7Z xzxz )i from i = rij + 2 to i = — 1 
(i.e. we skip the first exterior node and stop one node in from the outer 
boundary). We then return to the node we skipped over (i.e. i = rij + 1) and 
use cubic interpolation (using the the nodes nj, nj + 2, rij + 3 and rij + 4) to 
estimate (7Z XZXZ ) at that node. The Bianchi identity can then be applied to 
all the exterior nodes (except the node on the outer boundary). Finally, we 
use cubic extrapolation to compute the curvatures on the boundary nodes. 
This completes the computation of the curvatures. 



10.2 Maximal slicing 



The discrete form of the maximal lapse equation (8.1) is of the form 



= ai{N) i+1 + bi(N)i + Ci {N)i^ (10.5) 



for some set of coefficients a*, 5; and q (see Appendix [B] for the details). We 
wish to solve this set of equations subject to the following conditions 

= N >z at z = (10.6) 

= [N] at z = zj (10.7) 

= [N z ] at z = zj (10.8) 

1 = lim N (10.9) 

By reflection symmetry at z = we can easily extend the lattice to z < 0. 



Thus a discrete version of (10.6) would be iV_i = iV +1 . Continuity at z = zj 
allows us to use one value of N at zj, which we denote by Nj. However, the 
continuity of N tZ is not something we can prescribe but must be obtained by 
an iterative process (to be described below). We denote the left and right 
hand limits for iV z at z = zj by (iV 2 ) _ and (iV^) + respectively. We compute 
these one-sided limits, for a given set of (N)i, by a cubic extrapolation of 
(N jZ )i. This too requires some explanation. We start with the four nodes 
nearest to but excluding the junction. We then use cubic extrapolation of 
the (N)i on these nodes to extend the (N)i to the junction and two nodes 
beyond (we store this generated data in a separate array so as not to overwrite 
the data already defined on those nodes). Finally we use the standard non- 
uniform second order finite differences to estimate the first derivative at the 
junction. This computation is done twice, once for each one-sided limit. For 
the outer boundary we simply set 1 = iVoc,. 
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The discrete equations for (N)i are solved in three iterations with each iter- 
ation involving separate solutions for (iV)j on each side of the junction. The 
algorithm requires two guesses for N, one for iV at z = and one for N 
at z = zj which we denote by (GN)q and (GN)j respectively. With given 
values for these guesses we use a Thomas algorithm to solve the tri-diagonal 



system (10.5) for < z < zj and again for zj < z < z^. Our guesses 
are unlikely to be correct (at first) so we record the errors in the boundary 
conditions by E = (N) +1 - (JV)_i and Ej = (iV 2 ) + - (N z )~. Our aim is 
to choose the two guesses so that = E and = Ej. We chose three pairs 
of guesses (0, 0), (0, 1/2) and (1, 1) for ((GN) , (GN)j) and we recorded the 
corresponding errors as E { Q j) and E ( f for j = 1, 2, 3. Since the discrete equa- 
tions are linear and homogeneous in (N)i we can form a linear combination 
such as 

Ni = ai N\ l) + a 2 ivf } + a 3 N^ (10.10) 

to satisfy the boundary and junction by an appropriate choice of constants 
ax, «2 and a 3 . The result is a 3 by 3 system of equations 

1 = A^oo = oti + «2 + a 3 

= E = axE^ + a 2 ^ 2) + a 3 ^ 3) 

= Ej = ai Ej l) + a 2 Ef + a 3 Ej 3) 

which is easily solved for the three weights oti which in turn allows the final 



(correct) solution for the maximal lapse to be computed from (10.10). 



10.3 The time derivatives 

Spatial derivatives are calculated at each node using data from the surround- 
ing nodes, and in cases where this might draw in data from across the junc- 
tion, we first use cubic extrapolation to extend the data across the junction 
(which we store separately so as not to overwrite exiting data). 

With the exception of the junction node there is no ambiguity in applying the 
evolution equations to the nodes of the lattice. However, the discontinuities 
at the junction demand, once again, that we tread carefully near and at the 
junction. Consider JC ZZ which in geodesic slicing will be multiple valued at 
the junction. How do we handle this situation? We have already exhausted 



our supply of junction conditions in forming the two jump conditions (10.3) 



and (10.4) for the Riemann curvatures. So in the absence of any further 



information about K, xx we have no choice but to consider its left and right 
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hand limits as independent of each other (despite the loose coupling afforded 
by the evolution equations). Each term could be evolved by evaluating time 
derivatives built from one-sided limits of the source terms. There is however 
an easier approach which we found to work quite well. The idea is to re- 
interpret the junction node not as node on which to apply the evolution 
equations but rather as a convenient staging post to impose the junction 
conditions. In this view we do not evolve the data on the junction node. 
Rather we treat that data as kinematical which we compute by one-sided 
extrapolations of the surrounding data (which are evolved via the normal 
evolution equations). 



So in our code we use (4.3), (4.5), (4.6), (6.2) and (6.3) (subject to a minor 
change noted below) to evolve L xx , JC XX , 1C ZZ , v n and v z on the nodes i 
0,1,2, •• -nj—1 and i — nj+1, nj+2, rij+3, ■ ■ 



v n and v z 

— We use ( 4.4[ ) (again, 



see below) to evolve the L zz for all legs not connected to the junction. For the 
two legs attached to junction we use one-sided cubic extrapolation of dL zz /dt 
to compute their time derivatives. At the outer boundary we impose static 
boundary conditions for all of the data. 

There is one exception to this simple algorithm. We use a one-sided extrap- 
olation to set dL xx /dt at the node nj — 1. This proved to be essential for 
long term stability with maximal slicing (but made no difference in geodesic 
slicing). We can offer no reasonable explanation as to why this works other 
than the following admittedly vague rationalisation. By extrapolating the 
time derivatives outwards from the interior of the dust ball to the node nj — 1 
we might be halting or minimising the inward propagation of any errors that 
arise at the junction. Delving deeper into this mystery is best left for another 
time. 

There is one remaining subtlety that we must address. The careful reader 
may have noticed that in the present context we are treating the K, xx and 
}C ZZ as being defined on the nodes whereas the extrinsic curvatures arose in 
section [4] by approximating the integrals by a mid-point rule. Thus if we wish 
to use node based values for K xx and K zz we should use a Trapezoidal rule 
to estimate the integrals. This is a minor change and leads to the following 
node-based equations 



dL, 



dt 

dL z: 
dt 



- (NJC XX ) L xx + 

- {NK ZZ ) L zz + 



Nv 7 \ dL 



Nv z 



dz 



(10.11) 
(10.12) 



where the angle-brackets denotes an average of that quantity over the leg 
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while the square-brackets continues to denote the change across a leg. In 
fact for the L xx equation the angle-brackets are redundant (the end points 
carry identical values) but were retained simply for emphasis. Since the 
Riemann curvatures are already node-based we see that no such averaging 
is required for the extrinsic curvature equations (10.1) and (10.2). Note also 
that the spatial derivatives are also node based (by suitable choice of the 
finite difference operators). 



10.4 The spatial derivatives 



The evolution equations (4.3-4.6) and the momentum constraint (5.2) re- 
quire spatial derivatives of the L xx , L zz , K xx and K zz . For all but the two 
nodes either side of the junction (i.e. at nodes nj — 1 and nj + 1), and the 
junction itself, we employ second order non-uniform spatial derivatives as 
described in |Paper 1| On the two nodes either side of the junction we use 



one-sided quadratic extrapolation. This is the only point in the code where 
we used quadratic approximations and we do so because both linear and (in- 
terestingly) cubic interpolation lead to instabilities forming at the junction 
(at around t ~ 13 for cubic extrapolation and only for one of our models with 
nj = 240 and n M = 1200). The derivatives at the junction are computed 
last using one-sided cubic extrapolation. 

The only other spatial derivatives that need to be computed are the first 
and second derivatives of the lapse function (for use in the maximal slicing 



equation (8.1) and in the particle equations ( 6.2|6.3 )). Once again we use 
the second order non-uniform spatial derivatives from |Paper l| for all of the 



nodes with the exception of the five nodes centred on the junction. For nodes 
nj ± 2 and n j ± 1 we use cubic extrapolation to build an extended set data. 
This introduces some temporary and artificial nodes which we chose to be 
symmetric to the real nodes (e.g. when extending the data for node nj — 1 we 
create new nodes nj, nj + 1 that are the mirror images (in nj — 1) of nj — 2 
and nj — 3). The derivatives on nodes nj ± 2 and nj ± 1 are then computed 
on this extended data set using the standard non-uniform centred differences 
while the derivatives on the junction are computed using one-sided cubic 
extrapolation. 

Once the maximal slicing equation has been solved we do have the option of 
using that equation as an alternative way to calculate the second derivatives 
of the lapse. We chose not to do so because we did not want to give the 
smooth lattice method a helping hand - we want to test the method under 
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conditions closer (albeit in 1+1 form) to what we would expect for other 
spacetimes (i.e. for a true 3+1 evolution). 



10.5 The initial data 

We require two things of our initial data, first they must satisfy the con- 



straints (5.1) and (5.2), and second they must describe a time-symmetric 
initial slice. This last condition is readily satisfied upon setting K xx = 0, 
A^ 22 = and v n — 1, v z — which in turn ensures that the momentum con- 
straint is also satisfied. What we are left with is the Hamiltonian constraint, 
the leg lengths, L xx , L zz and the density p, i.e. we have one constraint for 
three (sets) of data. Clearly there are a range of options here, so what should 



we do? We turn once again to the scheme developed in Paper 1 There we 
chose to set the L zz and then use the Hamiltonian constraint to set the L xx . 
But here we also need the density. 

Keep in mind that our aim is neither to discover nor explore the Oppenheimer- 
Snyder solution but rather to use it as a test of the smooth lattice method. 
Thus it is not unreasonable to borrow some information from the exact solu- 
tion to set some of the data on the lattice, in particular the density. We recall 
here some basic equations from the exact solutions for the Oppenheimer- 
Snyder spacetime (see ED EH E2] ) • 

There are two free parameters in the solution, the ADM mass m and the 
Schwarzschild areal radius Rq of the dust ball. From these we can compute 
the proper radius of the dust ball zj, the FRW parameters a m and Xo > 
and the density p using 

o 2m 

sin 2 Xo = -- (10.13) 



2m 
sin 3 xo 



;i0.14) 



z.j = a m xo (10.15) 
8nkp = 4" (10.16) 

Clearly, we also have p = in the Schwarzschild exterior. 

We used these equations to set zj and p, for a given m and R . To this we 
added choices for the total length of the lattice z^, the number of interior 
nodes nj and the total number of nodes on the lattice. 
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Note that we still have the freedom to distribute the nodes along the z- 
axis (this amounts to setting the (L zz )i). We know that some of the spatial 
gradients are zero at z — 0, that they rise to a maximum near the junction 
and then settle down in the distant asymptotically flat regions of the lattice. 
Thus it makes sense to concentrate the nodes around the junction. With this 
in mind we chose to start at the junction and use a geometric progression to 
set the (L zz )i in both the interior and exterior regions. We chose the same 
geometric ratio in both regions while also requiring (L zz ) nj _i = (L zz ) nj . 
From here it is simple matter to compute all of the [L zz )i across the lattice. 

We now turn to the problem of setting L xx and the Riemann curvatures. 
By reworking the Hamiltonian constraint, geodesic deviation and Bianchi 
identity we find that across the lattice 



(10.17) 



(L X x)i — (L X x)i-l + JJ \ ((L xx )i-\ — (L xx )i_2) 

— ~(L zz )i_i {{L zz )i^i + [L zz )i-'i) (L xx TZ xzxz ) i l 

while the curvatures in the dust-ball are constant and are given by 

. . . . 8irkp 

(J^-xzxz)i (T^xyxyji g (10.18) 

and finally, in the Schwarzschild region, we find 

/ <3 ( K L XX ^^ ^ (.^ J xx^)i^\ 

\J^xzxz)i \J^xzxz!i—\ I r ( j 2 \ I T 2 ^ ) (10.19) 

\ L xx)i ~ \ L xx)i-\J 

{T^-xyxy)i ^{J^-xzxz) i (10.20) 

These equations can be used to set the L xx , lZ xyxy and lZ XZX z across the lattice 
(a process that will require the junction conditions for the curvatures). But 
to start the ball rolling we must make some choice for (L xx ) , and (L xx )i. 
Clearly (L xx )q = but for (L xz )i we are f ree to make any choice we like (we 
chose (L xx )i = 0.001(L 22 )o so that dL zz /dz = 0.001 at z — 0, as discussed 
below in section [Till). 



10.6 Density 

In section [7] we noted that the density can be computed using either the 



Hamiltonian constraint, in the form (7.1), or by the conservation equation 



(7.2). We find that, for long term stability when using the second method, 
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we are forced to use the Hamiltonian constraint at exactly the two nodes just 
inside the junction (i.e. at nodes n j — 2 and rij — 1). This was found by pure 
numerical experimentation. Why this should be so is unclear to us but it is 
probably tied to the same mechanism noted above (with regard to halting 
the inward propagation of errors from the junction by imposing "correct" 
values near the junction). 



11 Diagnostics 

From the known solution for the Oppenheimer- Snyder spacetime a number 
of useful diagnostics can be drawn. Here we will discuss those diagnostics 
which, in the following section, we will apply to our numerical results. 



For geodesic slicing it is rather easy to show [21] that the proper radius of 
the dust-ball zj varies with proper time t according to 

zj(t) = ^(l + cos V (t)) (11.1) 

where a m = 2m/ sin 3 xo and rj(t) is the solution of = —2t + a m (r] + sin?]) 
with 1] > (notice that Petrich et al. use rj where we use rj — tt). 

Another simple diagnostics arises from the central density which is given by 

p(t) = 24a 2 m (l + cosr ] (t))- 3 (11.2) 

This is singular when rj = n at which point the proper radius is zero and the 
dust ball has collapsed onto the singularity. This will occur after a proper 
time of 



and at this moment, or a short time before, we expect our code to crash. 

As the dust-ball collapses an outer apparent horizon will form and this too 
provides useful checks on our numerics. It is known that when the outer 
most apparent horizon forms it does so at the surface of the dust-ball. In 
our numerical code we locate the horizon by noting where on the radial 
axis the quantity dL xx /dz — K, XX L XX vanishes. The root of this equation is 
the location of the apparent horizon (this follows from the condition that 
= dA/du where A is the area of a 2-sphere and d/du is the outward 



pointing null vector to the 2-sphere, see Paper 1 for more details). The time 



at which the horizon forms is also well known and this affords yet another 
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check on our numerical results. For geodesic slicing it can be shown that the 
time, tg, and location zjf , of the apparent horizon are given by 

in 

-^(7T-2xo + sin(2 X o)) (11.4) 
sm d xo 

Tfl 

(l-cos(2 Xo )) (11.5) 

Note that in geodesic slicing the nodes are at rest relative to the Cauchy sur- 
faces and thus this time tg equals the proper time measured by the observer 
following that junction as it falls inwards and eventually meets the outward 
expanding event horizon. The quantity Zg measures the proper distance out 
from the centre of the dust-ball to the junction. 

Hawking's area theorem can also be used as a diagnostic. The theorem 
requires that the area of the event horizon should be constant once all of the 
dust has fallen within the event horizon. For our lattice this would require 
that the L xx on the event horizon should be constant for the remainder of 
the evolution. This is easily checked (by interpolating the values of L xx from 
the nodes onto the event horizon). 

Equations for the time and location of the horizon, as well as the density 
and radius diagnostics, are also available for maximal slicing but with one 
drawback - the equations as given by Petrich et al. require a numerical 
integration of some elliptic integrals. This introduces its own set of numerical 
issues and we found that our implementation of the Petrich equations could 
only be reliably used for t < 32 (for m = 1 and R — 5). Even so, this was 
sufficient time to allow for a useful comparison to be made. 

We also have one extra diagnostic for the case of maximal slicing. There it 
is known that the lapse function will, after an initial period, settle into an 
exponential decay. Petrich et al. show that N(t,0) ~ Aexp((3t) where A is 
a constant and f3 = -(2/3) (3/2) w -0.5443311. We can use this to test our 
code by measuring the slope of the logiV versus t. 

There are of course two other diagnostics - the Hamiltonian and momentum 
constraints. 

In summary we have the following set of diagnostics. 

• The constraints. 

• The history of the junction. 

• The history of the central density. 
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• The crash time for geodesic slicing. 

• The Petrich solution for maximal slicing. 

• The exponential collapse of the central lapse. 

• The time and location of the first apparent horizon. 

• The constancy of the area of the event horizon in the vacuum region. 

Clearly we have a raft of diagnostics and it is now time to turn to the actual 
results. 



12 Results 



Our aim was to write a code that used as few assumptions as needed to obtain 
reliable results. In the end we have split the computation of the lapse from 
the rest of the code. The evolution of the code takes as input (at each time 
step) the values of the lapse across the lattice. We do not use the Hamiltonian 
or momentum constraints apart from the two exceptions noted in sections 



10. 1| and |10.6[ We employ no artificial smoothing such as artificial viscosity 
nor do we add on any constraint preserving terms. Our time integrations 
are conducted using a 4th-order Runge-Kutta routine and our time step was 
updated after every time step by setting it equal to 1/2 the shortest L zz on 
the grid (which usually is the leg on or just inside the junction). This choice 
sets the Courant factor to 1/2 for legs near the junction (with smaller values 
for legs away from the junction). 

We set our initial data using S7ik — 1, m — 1, R = 5, Zoo = 400 and 
dL xx /dz = 0.001 at z — 0. We ran the code for three separate models, with 
(nj,noo) = (60,300), (120,600) and (240,1200) for both the geodesic and 
maximal slicing and one further model with (nj,noo) = (240,2400) for max- 
imal slicing. The results for a selection of quantities are displayed in Figures 
([5 -28). The first point to note is that the results are well behaved with no 



apparent instabilities even through to very late in the evolution. The junc- 
tion remains sharp without any noticeable smoothing and the constraints, 
though not zero, do not show the exponential growth often associated with 
unstable evolutions. 

We ran the geodesic code until it crashed at time t s = 12.41793 which 
compares well with the exact time = 12.41824 (note that the time step at 
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the crash was 8.19 x 10 6 which is considerably smaller than the initial time 
step of 5.38 x 1(T 3 ). 

For geodesic slicing we found the apparent horizon formed at t H = 10.87837 
and z H = 2.16534 while the exact values are = 10.87804 and z^ = 
2.1Qh27. While for maximal slicing the numerical values were t = 16.98238, 
z H = 2.38015 compared with the exact values tg = 17.02246, z% = 2.37971. 

For maximal slicing and the collapse of the lapse diagnostic we estimated the 
slope over the interval 25.0 < t < 35.0 and obtained (3 = —0.54424 compared 
with the exact value of —0.54433. 



In Figures ( 25p6 ) we have plotted the fractional errors in the radius and the 



central density for the first three models (as described above). For geodesic 
slicing the errors are very small. For maximal slicing the errors do decrease 
with increasing number of nodes however it would appear that the errors 
are not converging to zero. The simple explanation is that we set iV = 1 
on a finite outer boundary and this clearly incurs an error. To test this we 
re-ran our code with different choices for the location of the outer boundary 
(while retaining the same number of nodes). This showed that the peaks in 



Figures (26) varied inversely with the distance to the outer boundary z^. 
Incidentally, the broad peaks in those figures correspond to the formation of 
the apparent horizon. 

For maximal slicing we have taken a snapshot of the numerical data at a 



fixed time, Figure (24), to compare the density and the lapse with their 
exact values (from the Petrich code) across the lattice. Once again we see an 
initial convergence from coarse to fine resolutions but then the convergence 
appears to falter. This is also due to the use a of finite outer boundary, the 
peaks in the errors being proportional to Similar considerations apply 

to the snapshots of the Hamiltonian and momentum constraints, see Figure 



(23). The corresponding snapshot for geodesic slicing is shown in Figure (22). 
In this case the errors are not limited by z^ but instead depend only on rij 
and rioo and with the limited data available (only three models) it appears 
that the peaks in these figures reduce by a factor of about 4 for each doubling 

of (njjTloo). 



The fractional changes in the horizon L xx are shown in Figure (28). This 
shows that for t < 32 the horizon area varied by no more than 5 x 10 -2 
percent for the coarsest model improving to less than 1 x 10 -3 percent for 
the finest model. By t ~ 500 the error had grown to less than 2 percent for 
the finest model. 

We also ran our code using the Hamiltonian constraint to set the density and 
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found results very similar to those just given. 



13 Discussion 

The results just presented are very encouraging. They are consistent with 
our previous investigations of the smooth lattice method El E] yielding 
excellent results with only minor demands on computational resources. This 
gives us confidence that the method is viable but further tests are certainly 
required in particular an example in full 3+1 dimensions, without symme- 
tries, is imperative. This is a work in progress and we hope to report on this 
soon. 

One striking feature of the results for maximal slicing which we have so far 
ignored is the wave-like behaviour displayed in many of the plots (and similar 



behaviour was also noted in Paper 1). This is certainly not a gravitational 



wave (the spacetime is spherically symmetric). Can this behaviour be un- 
derstood from the evolution equations? Without delving too far into the 



analysis we note that the first order equations (4.3) and (4.5) can be recast 
as a single second order equation for L xx . This will involve d 2 L xx /dt 2 and 
the Riemann curvatures. But in these late times, where the waves are ap- 
parent, we see that \TZ XZXZ \ <§C [R- xyxy \ and thus the curvatures are dominated 
by 7Z xyxy which, through the geodesic deviation equation, (3.1), introduces 



d 2 L xx /dz 2 into the second order evolution equation for L xx . Thus we have in 
the one equation the two key elements of the one- dimensional wave equation 
for L xx and so wave-like behaviour is not surprising. Of course this is a very 
loose argument and there are many more terms to contend with before it can 
be said that the wave-like behaviour can be understood in standard terms. 
We will pursue this matter in a later paper. 



A The Darmois-Israel junction conditions 

Consider a spacetime (g,Ai) and let S be some 3-dimensional time like sur- 
face in M.. This surface will divide M. into two parts; one part, A4 C , to the 
left of S and another part, M n , to the right. In the absence of surface layers 
(e.g. infinitesimally thin shells of dust with non-zero energy) the Darmois- 
Israel junction conditions [T9l 120] ensure that g is a solution of Einstein's 
equations everywhere in M. provided it is a solution in Ai/S, and most im- 
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portantly, that the first and second fundamental forms on S are continuous 
across S. 



Suppose we denote the first and second fundamental forms on S by h and 
K respectively. Then each of these quantities can be calculated from the 
embedding of S in either Ai c or in J\A n . The junction conditions requires 
that both computations yield identical results, that is = [h] and = [K]. 

In our case we take (g,Ai) to be the Oppenheimer-Snyder spacetime and S 
to be the surface generated by the evolution of the surface of the dust. We 
will use a ~ symbol to denote quantities that live on S, for example, h and 
K will represent the 3-metric and extrinsic curvatures respectively on S. We 
extend this notation slightly to allow n to be unit (space like) normal to S 
in M. 

Our first task will be to express the junction conditions in terms of data on 
S. 

We know that L xx lies in S and thus the junction condition = [h] requires 
both = [L xx ] and = [dL xx /dt] while = [K] requires = [dL xx /dz] (note 
that d/dz is not normal to S but it can be resolved into pieces parallel and 
normal to S and the result follows). Looking back at the evolution equation 



(4.3) we see that this series of observations leads to the simple condition that 
= [ICxx] ■ We will make use of this result in the following discussions on the 
Riemann curvatures. Consider the Gauss equation for S, namely, 

where _L is the projection operator for S i.e. -L^ = 5 fJl u — n^h v . Since the 
vectors m^, are both tangent to S and since = [K^] we have 

= [ 4 R^ aul 3m^m x mym^} 

We can apply the Gauss equation once again, but this time for E rather than 
S, that is 

-L ( 4 Rfj,au/3) — R^iau/3 + K^ v K a p — K^K au 

This leads to the simple equation 

= [K xyxy ] (A.l) 

where we have used = [K xx ] and the fact that K^ v is diagonal. This is 
one of our two junction conditions for the Riemann curvature. The second 
condition will apply to 7Z XZXZ and as we shall soon see amounts to no more 
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than requiring continuity of the Hamiltonian constraint across the junction 
(as we would expect). 

We repeat the above procedure this time using the vectors and f 1 and 
after the first Gauss equation we find 

= fR^m^et^ 

Now t M is spanned by and m^, that is t M = v n n^ + v z m%, and thus we have 

= [vl± (R^pnW) m>f + 2v n v z ± (R^n") m u z m^ x 
+ v 2 z -\- (R^ap) m^m v z m a x mP x \ 

where we have also included the projection operator _L for £ (since and 
are both tangent to E) in preparation for the second application of the Gauss 
equation. This time we will need the Gauss equation and its contractions 
with n M , that is 

-L ( 4 R^au^) = R/j.au/3 + K^ u K a p — K^K^ 

-L ( 4 Rfj.au/3n 11 ) = K a fi\ v — K au \p 
JL ( 4 R^f3nW) = -± ( A R af3 ) + R aP + KK af3 - K ail K% 

Using the Einstein equations, 4 R af 3 = 87rk(T a p — (l/2)g a pT), the constraint 
equation — K^ u \ u = 87rk±(T^ u n u ) and the diagonal character of we 
find that 

-L (R^uap) m^ z m v z m a x mP x = K xzxz + K. XX JC ZZ 
-L (Rnvapn^) m v z m a x m? x = -4nkpv n v m 

-L {Rfi.avpn^n 1 ') m a z mP z = -Airkp + TZ xyxy + TZ XZXZ + JC XX K, ZZ + K? x 



XX 



and thus our junction condition can be reduced to 

1 







2^P^n T^xzxz K^xx^ 



(A.2) 



where we have used v\ — — 1 to eliminate v z . Looking back at our con- 
straint equations (5.1) we see that this last equation, along with = [TZ xyxy ] 
and = [/Ca-x], shows that the Hamiltonian constraint must be conserved 
across the junction (as expected). 
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B The maximal lapse equation 



Let (N)i be the node values of the lapse function across the lattice. Then 
using second order accurate finite differences (on a non-uniform grid) we 
obtain the following discrete equations 



with 



= a t (N) l+1 + 6,(7V). + Q (iV),-i (B.l) 



a t = [ + 1 ) (B.2) 



(L zz )i V. (L xx )i \ dz 
-2 / (L zz )i (dL, 



hi 



(L zz )i-i \(L xx )i \ dz J i 
4(L M ). ({^L zz ) t fdL 



1 (B.3) 



for z > and 



{L zz )i-i{L zz )i \ (L xx )i \ dz 
-2(L zz ).(R + nkp(8v 2 n + A)) (B.4) 



(B.5) 



(L zz )i (L zz )i_i 
-<L zz ) i 



(L zz )i-i(L zz )i 

for z — 0. 



2(Z^).(i? + 7rM8^ + 4)) (B.6) 



In the above equations we have introduced 2(L zz ) i = [L zz )i + {L Z z)%-\ and 
(A.L zz )i = {L zz )i — {L zz )i-\. 
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z 



Figure 1: In this figure we show how the lattice is constructed from two radial 
geodesies Go and G2 and the series of interconnecting legs L xx . The third geodesic 
Gi lies midway between Go and G2 is used to define the radial legs L zz . The grey 
patch to the left represents (part) of the dust ball. 




Figure 2: Here we display the 3-dimensional cell which we use to compute the 



energy density from the conservation equation ( 7.2 ) 
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t»6t = v^St 



= 



-v^n^r = N6t 



Figure 3: In this figure we introduce most of the kinematical quantities on the 
lattice. The dust particle's unit 4- velocity is v^, while 7^ is the drift vector, 
the unit normal to the Cauchy surface Ej, N is the lapse function and St is the 
proper time measured along the dust particle's trajectory. 
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Figure 4: In this diagram it is easy to see that [m^^] = dL xx /dz where 7 is the 
unit vector parallel to 7. This result is used in section [4] when deriving equation 
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Figure 5: The L xx and L zz leg lengths plotted from t = to t = 12 in steps of 0.8. 
The small dots denote the lattice node points. The larger diamonds denote the 
location of the apparent horizon. This occurs late in the evolution and appears 
only on the last two curves. The inward motion of the junction is also clearly 
evident in this plot. 
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Figure 6: The Riemann curvatures, lZ xyxy top and 1Z XZXZ bottom. Notice the 
flat profiles inside the dust ball. This feature can be seen in many of the following 
figures. 
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Figure 7: The extrinsic curvatures, fC xx top and IC ZZ bottom. 
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Figure 8: The constraints. These grow rapidly as the singularity is approached 
and this causes the first 10 curves to be too small to be seen on this scale. 
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Figure 9: The L xx leg lengths for t = to t = 32 in 20 steps (top, with < z < 10) 
and t = to t = 500 in 10 steps (bottom, with 10 < z < 400). This time 
the motion of the apparent horizon is much more noticeable than for the case of 
geodesic slicing. 
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Figure 10: The L zz leg lengths for t = to t = 32 in 20 steps (top, with 
< z < 10) and t = to t = 500 in 10 steps (bottom, with 10 < z < 400). Notice 
the extreme change in L zz at the junction. The curves bunch together late in the 



evolution due to the exponential collapse of the lapse (see Figure 17). 
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Figure 11: The Riemann curvature TZ xvxy . 
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Figure 12: The Riemann curvature TZ XZXZ . 
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Figure 13: The extrinsic curvature K, xx . 
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Figure 14: The extrinsic curvature 1C ZZ . 
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Figure 15: The Hamiltonian constraint. This shows a slowly growing error for 
later times. 
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Figure 16: The momentum constraint. This shows a similar slow growing peak 
as seen in the Hamiltonian constraint. 
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Figure 17: A logarithmic plot of the lapse. The even gaps between the curves 
shows clearly that the collapse is exponential in time. Note the extreme value of 
the lapse at the origin for late times, of order 10~ 110 . 
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Figure 18: The trace of the extrinsic curvature K. This should be zero for all 
time. In the lower plot the lapse has collapsed in that part of the lattice and thus 
there is no apparent evolution in K. In the upper plot {nj^n^) = (240,1200) 
while for lower plot we used (nj,noo) = (240,2400). 
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Figure 19: The particle velocity v. 
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Figure 20: The particle velocity v z . 
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Figure 21: The The rest density for geodesic slicing (top) and maximal slicing 
(bottom). 
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Figure 22: This is a snapshot of the constraints across the lattice at a fixed time 
in geodesic slicing. The three curves correspond to the three models described in 
the text. Note that the peaks decrease rapidly as the number of lattice nodes is 
increased. The horizontal axes have been truncated to give a better view of the 
data. 
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Figure 23: The Hamiltonian and momentum constraints across the lattice at a 
fixed time with maximal slicing. The peak occurs near the junction and appears 
to vary as l/n^. 
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Figure 24: The fractional errors for the lapse and the rest density across the 
lattice with maximal slicing. The large error in the coarsest lattice is probably- 
due to having two few nodes. The finer lattice show much better errors but note 
that the lapse appears not to converge at the origin. This is due to the use of a 
finite outer boundary for the lapse. 



54 



5.00e-08 



0.00e+00 



-5.00e-08 - 



-1.00e-07 - 



-1.50e-07 



g -2.00e-07 
O 



-2.50e-07 - 



-3.00e-07 




0.0 



2.0 



4.0 6.0 8.0 10.0 

Proper time r 



12.0 



14.0 




-2.00e-10 



4.0 6.0 8.0 10.0 

Proper time r 



12.0 14.0 



Figure 25: The fractional errors for the radius and the central density for geodesic 
slicing. The convergence is clear and it is rapid (we make no attempt to estimate 
the order of the convergence) . 
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Figure 26: The fractional errors for the radius and the central density for maximal 
slicing. The peak occurs around the time when the apparent horizon forms. The 
height of the peak for the finest resolution is limited by the location of the outer 
boundary. Doubling halves the height of the peak. 
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Figure 27: The lapse at the origin for the three models superimposed on the 
exact data of Petrich et al. (top) and for the single long term model (bottom). 
This shows clearly that the lapse collapses exponentially. 



57 




16.0 18.0 20.0 22.0 24.0 26.0 28.0 30.0 32.0 34.0 
Proper time r 




50.0 100.0 150.0 200.0 250.0 300.0 350.0 400.0 450.0 500.0 
Proper time r 



Figure 28: The fractional change in the L xx on the horizon for the three models 
(top) and for the long-term integration (bottom) in maximal slicing). These errors 
should be zero by Hawking area theorem. 
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